Masterthesis script KMI
Heleen Van Ransbeeck
2023-02-27
Data inladen
Baits<-read.table(file = "~/GitHub/vespa_analyses/Input/Baits_useful.txt", header=TRUE, sep="\t", na.strings=c(""," ","NA"))
Flights<-read.table(file = "~/GitHub/vespa_analyses/Input/Flights_KMI.txt", header=TRUE, sep="\t", na.strings=c(""," ","NA"))
Nests<-read.table(file = "~/GitHub/vespa_analyses/Input/Nests_useful.txt", header=TRUE, sep="\t", na.strings=c(""," ","NA"))
Shortind_incomplete<-read.table(file = "~/GitHub/vespa_analyses/Input/Shortind_incomplete.txt", header=TRUE, sep="\t", na.strings=c(""," ","NA"))
Shortpot_incomplete<-read.table(file = "~/GitHub/vespa_analyses/Input/Shortpot_incomplete.txt", header=TRUE, sep="\t", na.strings=c(""," ","NA"))
Individuals<-read.table(file = "~/GitHub/vespa_analyses/Input/Individuals.txt", header=TRUE, sep="\t", na.strings=c(""," ","NA"))Aanpassingen datasets
# Overbodige kolommen wegdoen, foutje in excel
Flights<-Flights[1:(length(Flights)-1)]
Shortind_incomplete<-Shortind_incomplete[1:(length(Shortind_incomplete)-2)]
Shortpot_incomplete<-Shortpot_incomplete[1:(length(Shortpot_incomplete)-2)]
#Foerageersnelheid en vliegsnelheid (m/s) toevoegen aan tabel
Flights$test_rule_of_thumb<-Flights$Distance/(Flights$Flighttime_min*60)
Shortind_incomplete$test_rule_of_thumb<-Shortind_incomplete$Distance/(Shortind_incomplete$Flighttime_min*60)
Shortpot_incomplete$test_rule_of_thumb<-Shortpot_incomplete$Distance/(Shortpot_incomplete$Flighttime_min*60)
Flights$ForagingSpeed<-Flights$Distance*2/(Flights$Flighttime_min*60)
Shortind_incomplete$ForagingSpeed<-Shortind_incomplete$Distance*2/(Shortind_incomplete$Flighttime_min*60)
Shortpot_incomplete$ForagingSpeed<-Shortpot_incomplete$Distance*2/(Shortpot_incomplete$Flighttime_min*60)
## Urbanisatiecirkels per vliegtijd toevoegen
Flights<-merge(Flights, Baits[, c("NestID", "BaitID", "Urbanisation25m", "Urbanisation50m", "Urbanisation100m")], by=c("NestID", "BaitID"), all=TRUE)
Shortind_incomplete<-merge(Shortind_incomplete, Baits[, c("NestID", "BaitID", "Urbanisation25m", "Urbanisation50m", "Urbanisation100m")], by=c("NestID", "BaitID"))
Shortpot_incomplete<-merge(Shortpot_incomplete, Baits[, c("NestID", "BaitID", "Urbanisation25m", "Urbanisation50m", "Urbanisation100m")], by=c("NestID", "BaitID"))
## Urbanisatietrajecten per vliegtijd toevoegen
Flights<-merge(Flights, Baits[, c("NestID", "BaitID", "Traject25m", "Traject50m", "Traject100m")], by=c("NestID", "BaitID"), all=TRUE)
Shortind_incomplete<-merge(Shortind_incomplete, Baits[, c("NestID", "BaitID", "Traject25m", "Traject50m", "Traject100m")], by=c("NestID", "BaitID"))
Shortpot_incomplete<-merge(Shortpot_incomplete, Baits[, c("NestID", "BaitID", "Traject25m", "Traject50m", "Traject100m")], by=c("NestID", "BaitID"))
## Gewicht individu per vliegtijd toevoegen
Flights<-merge(Flights, Individuals[, c("NestID", "BaitID","ColorInd", "Weight_ind")], by=c("NestID", "BaitID", "ColorInd"), all=TRUE)
Shortind_incomplete<-merge(Shortind_incomplete, Individuals[, c("NestID", "BaitID", "ColorInd", "Weight_ind")], by=c("NestID", "BaitID", "ColorInd"))
Shortpot_incomplete<-merge(Shortpot_incomplete, Individuals[, c("NestID", "BaitID", "ColorInd","Weight_ind")], by=c("NestID", "BaitID", "ColorInd"))
## Nesthoogte per vliegtijd toevoegen
Flights<-merge(Flights, Nests[, c("NestID", "Height")], by=c("NestID"), all=TRUE)
Shortind_incomplete<-merge(Shortind_incomplete, Nests[, c("NestID", "Height")], by=c("NestID"))
Shortpot_incomplete<-merge(Shortpot_incomplete, Nests[, c("NestID", "Height")], by=c("NestID"))
library(ggplot2)
library(ggpubr)Volgende datasets zijn gecreëerd:
- Flights: dataset met alle vliegtijden
- Shortind_incomplete: dataset met kortste vliegtijden per individu
- Shortpot_incomplete: dataset met kortste vliegtijden per wiekpot
Voor het model: Flight time ~ Distance gebruiken we de dataset vliegtijden per individu Omdat we hiermee de theoretische regel 1min=100m kunnen verifiëren. De imkers nemen hiervoor altijd kortste meting
Voor modellen met weerparameters, gewicht, urbanisatie: Hiervoor nemen we telkens de hele dataset, omdat elke meting van deze factoren afhangt.
Datapunten op 2km weglaten
We gaan verder met de dataset zonder de outliers in Melle
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionFlights2<-Flights %>% filter(NestID!=1)
Shortind2_incomplete<-Shortind_incomplete %>% filter(NestID!=1)
Shortpot2_incomplete<-Shortpot_incomplete %>% filter(NestID!=1)
# ID maken per individu
Flights2$Individualcode <- paste(Flights2$NestID, Flights2$BaitID,Flights2$ColorInd, sep="_")
Shortind2_incomplete$Individualcode <- paste(Shortind2_incomplete$NestID, Shortind2_incomplete$BaitID,Shortind2_incomplete$ColorInd, sep="_")VRAAG: Hoe grafiek maken van parameters mixed model?
1.1) ForagingSpeed vs Temperature
Graphs
Model Output
Normality niet ok!
Niet significant!
library(lmerTest)## Loading required package: lme4## Loading required package: Matrix##
## Attaching package: 'lmerTest'## The following object is masked from 'package:lme4':
##
## lmer## The following object is masked from 'package:stats':
##
## steplibrary(lme4)
model_temp_all<-lmer(ForagingSpeed ~ Temperature_KMI + (1|NestID) + (1|NestID:BaitID) + (1|Individualcode), na.action=na.omit, data=Flights2)
summary(model_temp_all)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula:
## ForagingSpeed ~ Temperature_KMI + (1 | NestID) + (1 | NestID:BaitID) +
## (1 | Individualcode)
## Data: Flights2
##
## REML criterion at convergence: 595.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.19526 -0.45657 0.04751 0.38435 3.07187
##
## Random effects:
## Groups Name Variance Std.Dev.
## Individualcode (Intercept) 0.1802 0.4245
## NestID:BaitID (Intercept) 0.7973 0.8929
## NestID (Intercept) 0.7893 0.8884
## Residual 0.3540 0.5949
## Number of obs: 230, groups: Individualcode, 91; NestID:BaitID, 68; NestID, 34
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 3.234075 0.537550 135.595552 6.016 1.57e-08 ***
## Temperature_KMI -0.001193 0.026258 175.559355 -0.045 0.964
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## Temprtr_KMI -0.919anova(model_temp_all, ddf="Satterthwaite", type=3)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Temperature_KMI 0.00073042 0.00073042 1 175.56 0.0021 0.9638res<-residuals(model_temp_all)
shapiro.test(res)##
## Shapiro-Wilk normality test
##
## data: res
## W = 0.98221, p-value = 0.005473qqnorm(res)
qqline(res)
1.2) Flight time vs Distance + Temperature
Graph
Klopt dit? Kan dit ook met mixed model?
library(knitr)
library(kableExtra)## Warning in !is.null(rmarkdown::metadata$output) && rmarkdown::metadata$output
## %in% : 'length(x) = 2 > 1' in coercion to 'logical(1)'##
## Attaching package: 'kableExtra'## The following object is masked from 'package:dplyr':
##
## group_rowslibrary("plot3D")
x <- Flights2$Distance
y <- Flights2$Temperature_KMI
z <- Flights2$Flighttime_min
fit <- lm(z ~ x + y, na.action=na.exclude)
x.pred <- seq(min(x[!is.na(x)]), max(x[!is.na(x)]), length.out = 20)
y.pred <- seq(min(y[!is.na(y)]), max(y[!is.na(y)]), length.out = 20)
xy <- expand.grid( x = x.pred, y = y.pred)
z.pred <- matrix(predict(fit, newdata = xy),
nrow = 20, ncol = 20)
fitpoints <- predict(fit)
scatter3D(x, y, z, pch = 19, cex = 0.6, colvar=FALSE, col="dodgerblue3", theta = 210, phi = 10, bty="u", col.panel ="grey93", expand =0.4, col.grid = "white", xlab = "Distance", ylab = "Temperature", zlab = "Flight time", surf = list(x = x.pred, y = y.pred, z = z.pred,
facets = TRUE, col=ramp.col(col = c("dodgerblue4", "seagreen2"), n = 100, alpha=0.8), fit = fitpoints, border="black"),main = "Flight time vs Distance + Temperature")
Met plotly
Had met deze link problemen met expand.grid https://stackoverflow.com/questions/38331198/add-regression-plane-to-3d-scatter-plot-in-plotly Dan maar grid dat van hierboven gebruikt, maar nu klopt er precies iets niet? Alle datapunten liggen boven het oppervlak.
library(plotly)##
## Attaching package: 'plotly'## The following object is masked from 'package:ggplot2':
##
## last_plot## The following object is masked from 'package:stats':
##
## filter## The following object is masked from 'package:graphics':
##
## layoutlibrary(tidyverse)## ── Attaching packages
## ───────────────────────────────────────
## tidyverse 1.3.2 ──## ✔ tibble 3.1.8 ✔ purrr 1.0.1
## ✔ tidyr 1.2.1 ✔ stringr 1.5.0
## ✔ readr 2.1.4 ✔ forcats 1.0.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ tidyr::expand() masks Matrix::expand()
## ✖ plotly::filter() masks dplyr::filter(), stats::filter()
## ✖ kableExtra::group_rows() masks dplyr::group_rows()
## ✖ dplyr::lag() masks stats::lag()
## ✖ tidyr::pack() masks Matrix::pack()
## ✖ tidyr::unpack() masks Matrix::unpack()FlightsNoNA<-Flights2 %>% filter(Flighttime_min!="NA")
FlightsNoNA<-FlightsNoNA %>% filter(Temperature_KMI!="NA")
fig <- plot_ly(FlightsNoNA, x = ~Distance, y = ~Temperature_KMI, z = ~Flighttime_min, size=1)
fig <- fig %>% add_markers()
fig <- fig %>% layout(scene = list(xaxis = list(title = 'Distance'),
yaxis = list(title = 'Temperature_KMI'),
zaxis = list(title = 'Flight time (min)')))
p2 <- add_trace(p = fig,
z = z.pred,
x = seq(2100, 0, by = -100),
y = seq(0, 30, by = 10),
type = "surface")
p2Model Ouput
Normality niet ok!
Niet significant
model_tempdist_all<-lmer(Flighttime_min ~ Distance + Temperature_KMI + (1|NestID) + (1|NestID:BaitID) + (1|Individualcode), na.action=na.exclude, data=Flights2)## boundary (singular) fit: see help('isSingular')summary(model_tempdist_all)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Flighttime_min ~ Distance + Temperature_KMI + (1 | NestID) +
## (1 | NestID:BaitID) + (1 | Individualcode)
## Data: Flights2
##
## REML criterion at convergence: 1005.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.2425 -0.4262 -0.1822 0.1635 4.4765
##
## Random effects:
## Groups Name Variance Std.Dev.
## Individualcode (Intercept) 0.1443 0.3799
## NestID:BaitID (Intercept) 1.8910 1.3751
## NestID (Intercept) 0.0000 0.0000
## Residual 3.2814 1.8115
## Number of obs: 230, groups: Individualcode, 91; NestID:BaitID, 68; NestID, 34
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 3.5399788 0.9074782 66.9650223 3.901 0.000225 ***
## Distance 0.0056150 0.0008062 62.3323406 6.965 2.42e-09 ***
## Temperature_KMI -0.0686672 0.0445457 91.2883267 -1.541 0.126655
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Distnc
## Distance -0.336
## Temprtr_KMI -0.903 -0.010
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')anova(model_tempdist_all, ddf="Satterthwaite", type=3)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Distance 159.173 159.173 1 62.332 48.5081 2.415e-09 ***
## Temperature_KMI 7.797 7.797 1 91.288 2.3762 0.1267
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1res<-residuals(model_tempdist_all)
shapiro.test(res)##
## Shapiro-Wilk normality test
##
## data: res
## W = 0.7916, p-value < 2.2e-16qqnorm(res)
qqline(res)
2) ForagingSpeed vs Cloudcoverage
Graphs
ggplot(Flights2, aes(x=Cloudcoverage_KMI, y=ForagingSpeed)) + geom_point(color="lightblue") + geom_smooth(method="lm", formula =y ~ x, se=TRUE, fullrange=FALSE, level=0.95, color="darkblue") + ggtitle("All data KMI")+ theme(plot.title = element_text(hjust = 0.5, size=10))
Model Ouput
Licht significant
Normality niet ok!
model_cloud_all<-lmer(ForagingSpeed ~ Cloudcoverage_KMI + (1|NestID) + (1|NestID:BaitID) + (1|Individualcode), na.action=na.omit, data=Flights2)
summary(model_cloud_all)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula:
## ForagingSpeed ~ Cloudcoverage_KMI + (1 | NestID) + (1 | NestID:BaitID) +
## (1 | Individualcode)
## Data: Flights2
##
## REML criterion at convergence: 491.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.2243 -0.5101 0.0644 0.4260 2.9113
##
## Random effects:
## Groups Name Variance Std.Dev.
## Individualcode (Intercept) 0.1866 0.4320
## NestID:BaitID (Intercept) 0.9826 0.9912
## NestID (Intercept) 0.7492 0.8655
## Residual 0.3547 0.5956
## Number of obs: 191, groups: Individualcode, 74; NestID:BaitID, 52; NestID, 23
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 3.2816 0.2715 24.9103 12.088 6.45e-12 ***
## Cloudcoverage_KMI -0.4237 0.2090 141.7036 -2.027 0.0445 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## Cldcvrg_KMI -0.326anova(model_cloud_all, ddf="Satterthwaite", type=3)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Cloudcoverage_KMI 1.4581 1.4581 1 141.7 4.1106 0.04449 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1res<-residuals(model_cloud_all)
shapiro.test(res)##
## Shapiro-Wilk normality test
##
## data: res
## W = 0.98299, p-value = 0.0204qqnorm(res)
qqline(res)
3.1) ForagingSpeed vs Windspeed
Graphs
Telkens voor de modellen:
ForagingSpeed ~ Windspeed
ForagingSpeed ~ Windspeed²
plot19<-ggplot(Flights2, aes(x=WindSpeed_KMI, y=ForagingSpeed)) + geom_point(color="khaki3") + geom_smooth(method="lm", formula =y ~ x, se=TRUE, fullrange=FALSE, level=0.95, color="olivedrab") + ggtitle("ForagingSpeed ~ Windspeed KMI")+ theme(plot.title = element_text(hjust = 0.5, size=10))
plot22<-ggplot(Flights2, aes(x=WindSpeed_KMI, y=ForagingSpeed)) + geom_point(color="khaki3") + geom_smooth(method="lm", formula =y ~ I(x^2), se=TRUE, fullrange=FALSE, level=0.95, color="olivedrab") + ggtitle("ForagingSpeed ~ Windspeed² KMI")+ theme(plot.title = element_text(hjust = 0.5, size=10))
ggarrange(plot19, plot22 + rremove("x.text"),
labels = c("A", "B"),
ncol = 2, nrow = 1)
Model Output
Beide modellen hoogsignificant
Normality voor beiden ok!
model_wind_all<-lmer(ForagingSpeed ~ WindSpeed_KMI + (1|NestID) + (1|NestID:BaitID) + (1|Individualcode), na.action=na.omit, data=Flights2)
summary(model_wind_all)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: ForagingSpeed ~ WindSpeed_KMI + (1 | NestID) + (1 | NestID:BaitID) +
## (1 | Individualcode)
## Data: Flights2
##
## REML criterion at convergence: 572.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.55366 -0.52807 0.07123 0.41517 2.99472
##
## Random effects:
## Groups Name Variance Std.Dev.
## Individualcode (Intercept) 0.1681 0.4100
## NestID:BaitID (Intercept) 0.7091 0.8421
## NestID (Intercept) 0.8242 0.9079
## Residual 0.3160 0.5621
## Number of obs: 230, groups: Individualcode, 91; NestID:BaitID, 68; NestID, 34
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 3.94765 0.25805 70.73569 15.30 < 2e-16 ***
## WindSpeed_KMI -0.19399 0.03984 165.66865 -4.87 2.59e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## WindSpd_KMI -0.584anova(model_wind_all, ddf="Satterthwaite", type=3)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## WindSpeed_KMI 7.4929 7.4929 1 165.67 23.713 2.594e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1model_wind2_all<-lmer(ForagingSpeed ~ I(WindSpeed_KMI^2) + (1|NestID) + (1|NestID:BaitID) + (1|Individualcode), na.action=na.omit, data=Flights2)
summary(model_wind2_all)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula:
## ForagingSpeed ~ I(WindSpeed_KMI^2) + (1 | NestID) + (1 | NestID:BaitID) +
## (1 | Individualcode)
## Data: Flights2
##
## REML criterion at convergence: 580
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.26566 -0.54124 0.08453 0.41515 2.93663
##
## Random effects:
## Groups Name Variance Std.Dev.
## Individualcode (Intercept) 0.1414 0.3760
## NestID:BaitID (Intercept) 0.7246 0.8513
## NestID (Intercept) 0.8080 0.8989
## Residual 0.3295 0.5740
## Number of obs: 230, groups: Individualcode, 91; NestID:BaitID, 68; NestID, 34
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 3.599231 0.226024 45.995965 15.92 < 2e-16 ***
## I(WindSpeed_KMI^2) -0.022824 0.005175 191.376210 -4.41 1.72e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## I(WS_KMI^2) -0.387anova(model_wind2_all, ddf="Satterthwaite", type=3)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## I(WindSpeed_KMI^2) 6.4084 6.4084 1 191.38 19.451 1.719e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1res<-residuals(model_wind_all)
shapiro.test(res)##
## Shapiro-Wilk normality test
##
## data: res
## W = 0.98923, p-value = 0.0834qqnorm(res)
qqline(res)
res<-residuals(model_wind2_all)
shapiro.test(res)##
## Shapiro-Wilk normality test
##
## data: res
## W = 0.98929, p-value = 0.08548qqnorm(res)
qqline(res)
3.2) ForagingSpeed vs Windspeed | afhankelijk van de windrichting
Hiervoor werd telkens voor elke meting nagegaan of de hoornaar met meewind (tailwind), tegenwind (upwind) of loodrechte wind (perpendicular) te maken had. Dit volgens de formules:
|𝜃flight −𝜃wind| ≤ 45 is tailwind
45 < |𝜃flight −𝜃wind|<135 is (quasi) perpendicular
|𝜃 flight −𝜃wind| ≥135 upwind
Graphs
ggplot(data=subset(Flights2, !is.na(Wind)), aes(x= Wind, col=Wind, y=ForagingSpeed)) + geom_boxplot()## Warning: Removed 14 rows containing non-finite values (`stat_boxplot()`).
ggplot(data=subset(Flights2, !is.na(Wind)), aes(x= WindSpeed_KMI, col=Wind, y=ForagingSpeed)) + geom_point() + facet_grid(~Wind) + geom_smooth(method="lm", formula = y ~ I(x^2), se=TRUE, fullrange=FALSE, level=0.95)## Warning: Removed 14 rows containing non-finite values (`stat_smooth()`).## Warning: Removed 14 rows containing missing values (`geom_point()`).
Model Output
Anova van Foragingspeed en de 3 windinvloeden
- interpretatie?
windanova<-lmer(ForagingSpeed ~ Wind + (1|NestID) + (1|NestID:BaitID) + (1|Individualcode), na.action=na.exclude, data=Flights2)
summary(windanova)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: ForagingSpeed ~ Wind + (1 | NestID) + (1 | NestID:BaitID) + (1 |
## Individualcode)
## Data: Flights2
##
## REML criterion at convergence: 592.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.1864 -0.4629 0.0327 0.3908 3.0729
##
## Random effects:
## Groups Name Variance Std.Dev.
## Individualcode (Intercept) 0.1877 0.4332
## NestID:BaitID (Intercept) 0.7814 0.8840
## NestID (Intercept) 0.7900 0.8888
## Residual 0.3548 0.5956
## Number of obs: 230, groups: Individualcode, 91; NestID:BaitID, 68; NestID, 34
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 3.24826 0.22447 40.83023 14.471 <2e-16 ***
## Windtailwind -0.14580 0.19936 214.92685 -0.731 0.465
## Windupwind -0.02152 0.18402 182.45501 -0.117 0.907
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Wndtlw
## Windtailwnd -0.240
## Windupwind -0.279 0.224ForagingSpeed vs WindSpeed voor Tailwind dataset
Hoogsignificant
Normality ok!
data_tailwind<-Flights2 %>% filter(Wind!="tailwind")
model_tailwind<-lmer(ForagingSpeed ~ WindSpeed_KMI + (1|NestID) + (1|NestID:BaitID) + (1|Individualcode), na.action=na.omit, data=data_tailwind)
summary(model_tailwind)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: ForagingSpeed ~ WindSpeed_KMI + (1 | NestID) + (1 | NestID:BaitID) +
## (1 | Individualcode)
## Data: data_tailwind
##
## REML criterion at convergence: 476.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.45714 -0.55810 0.04615 0.41914 2.91009
##
## Random effects:
## Groups Name Variance Std.Dev.
## Individualcode (Intercept) 0.1738 0.4169
## NestID:BaitID (Intercept) 0.7757 0.8807
## NestID (Intercept) 0.7954 0.8919
## Residual 0.3318 0.5760
## Number of obs: 188, groups: Individualcode, 75; NestID:BaitID, 54; NestID, 31
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 3.81145 0.27541 60.89085 13.839 < 2e-16 ***
## WindSpeed_KMI -0.16562 0.04216 130.28317 -3.929 0.000138 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## WindSpd_KMI -0.574anova(model_tailwind, ddf="Satterthwaite", type=3)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## WindSpeed_KMI 5.121 5.121 1 130.28 15.434 0.0001379 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1res<-residuals(model_tailwind)
shapiro.test(res)##
## Shapiro-Wilk normality test
##
## data: res
## W = 0.98883, p-value = 0.1474qqnorm(res)
qqline(res)
ForagingSpeed vs WindSpeed voor Perpendicular dataset
Niet significant
Normality niet ok!
data_perpendicular<-Flights2 %>% filter(Wind!="perpendicular")
model_perpendicular<-lmer(ForagingSpeed ~ WindSpeed_KMI + (1|NestID) + (1|NestID:BaitID) + (1|Individualcode), na.action=na.omit, data=data_perpendicular)
summary(model_perpendicular)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: ForagingSpeed ~ WindSpeed_KMI + (1 | NestID) + (1 | NestID:BaitID) +
## (1 | Individualcode)
## Data: data_perpendicular
##
## REML criterion at convergence: 259.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.4499 -0.3729 0.0493 0.4221 3.3570
##
## Random effects:
## Groups Name Variance Std.Dev.
## Individualcode (Intercept) 0.01943 0.1394
## NestID:BaitID (Intercept) 1.12489 1.0606
## NestID (Intercept) 0.86845 0.9319
## Residual 0.28222 0.5312
## Number of obs: 98, groups: Individualcode, 47; NestID:BaitID, 40; NestID, 27
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 3.41766 0.36128 50.36655 9.460 9.37e-13 ***
## WindSpeed_KMI -0.09641 0.06997 48.59291 -1.378 0.175
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## WindSpd_KMI -0.673anova(model_perpendicular, ddf="Satterthwaite", type=3)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## WindSpeed_KMI 0.53583 0.53583 1 48.593 1.8986 0.1745res<-residuals(model_perpendicular)
shapiro.test(res)##
## Shapiro-Wilk normality test
##
## data: res
## W = 0.95086, p-value = 0.001085qqnorm(res)
qqline(res)
ForagingSpeed vs WindSpeed voor Upwind dataset
Hoogsignificant
Normality ok!
data_upwind<-Flights2 %>% filter(Wind!="upwind")
model_upwind<-lmer(ForagingSpeed ~ WindSpeed_KMI + (1|NestID) + (1|NestID:BaitID) + (1|Individualcode), na.action=na.omit, data=data_upwind)
summary(model_upwind)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: ForagingSpeed ~ WindSpeed_KMI + (1 | NestID) + (1 | NestID:BaitID) +
## (1 | Individualcode)
## Data: data_upwind
##
## REML criterion at convergence: 432.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.6138 -0.4276 0.1099 0.4081 2.2258
##
## Random effects:
## Groups Name Variance Std.Dev.
## Individualcode (Intercept) 0.1331 0.3648
## NestID:BaitID (Intercept) 0.8066 0.8981
## NestID (Intercept) 0.5930 0.7700
## Residual 0.3105 0.5572
## Number of obs: 174, groups: Individualcode, 71; NestID:BaitID, 56; NestID, 28
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 4.13221 0.27666 64.34513 14.936 < 2e-16 ***
## WindSpeed_KMI -0.22679 0.04392 122.09885 -5.163 9.55e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## WindSpd_KMI -0.631anova(model_upwind, ddf="Satterthwaite", type=3)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## WindSpeed_KMI 8.2778 8.2778 1 122.1 26.659 9.555e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1res<-residuals(model_upwind)
shapiro.test(res)##
## Shapiro-Wilk normality test
##
## data: res
## W = 0.98133, p-value = 0.01951qqnorm(res)
qqline(res)
Model Selection
fullmodel<-lmer(ForagingSpeed ~ Traject100m + WindSpeed_KMI + Temperature_KMI + Cloudcoverage_KMI + Weight_ind + (1|NestID) + (1|NestID:BaitID) + (1|Individualcode), na.action=na.omit, data=Flights2)## boundary (singular) fit: see help('isSingular')summary(fullmodel)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: ForagingSpeed ~ Traject100m + WindSpeed_KMI + Temperature_KMI +
## Cloudcoverage_KMI + Weight_ind + (1 | NestID) + (1 | NestID:BaitID) +
## (1 | Individualcode)
## Data: Flights2
##
## REML criterion at convergence: 132.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.4288 -0.5430 0.0748 0.5100 1.6174
##
## Random effects:
## Groups Name Variance Std.Dev.
## Individualcode (Intercept) 0.07951 0.282
## NestID:BaitID (Intercept) 0.66749 0.817
## NestID (Intercept) 0.00000 0.000
## Residual 0.25708 0.507
## Number of obs: 71, groups: Individualcode, 16; NestID:BaitID, 12; NestID, 7
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 5.24098 2.24293 19.02657 2.337 0.0305 *
## Traject100m -5.85766 2.63126 8.90757 -2.226 0.0533 .
## WindSpeed_KMI -0.15017 0.20012 15.15305 -0.750 0.4645
## Temperature_KMI -0.05689 0.09633 21.97521 -0.591 0.5608
## Cloudcoverage_KMI 0.59392 0.44256 38.57798 1.342 0.1874
## Weight_ind -1.12665 3.95322 11.66851 -0.285 0.7806
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Trj100 WS_KMI Tm_KMI Cl_KMI
## Traject100m -0.495
## WindSpd_KMI 0.684 -0.366
## Temprtr_KMI -0.809 0.309 -0.857
## Cldcvrg_KMI -0.366 0.278 -0.308 0.397
## Weight_ind -0.794 0.146 -0.403 0.418 0.034
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')model1<-lmer(ForagingSpeed ~ Traject100m + WindSpeed_KMI + Temperature_KMI + Cloudcoverage_KMI + (1|NestID) + (1|NestID:BaitID) + (1|Individualcode), na.action=na.omit, data=Flights2)
summary(model1)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: ForagingSpeed ~ Traject100m + WindSpeed_KMI + Temperature_KMI +
## Cloudcoverage_KMI + (1 | NestID) + (1 | NestID:BaitID) +
## (1 | Individualcode)
## Data: Flights2
##
## REML criterion at convergence: 469.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.44419 -0.54207 0.03008 0.46341 2.91374
##
## Random effects:
## Groups Name Variance Std.Dev.
## Individualcode (Intercept) 0.1602 0.4002
## NestID:BaitID (Intercept) 0.9180 0.9581
## NestID (Intercept) 0.2815 0.5306
## Residual 0.3234 0.5687
## Number of obs: 191, groups: Individualcode, 74; NestID:BaitID, 52; NestID, 23
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 4.06254 0.74726 64.14701 5.437 9.02e-07 ***
## Traject100m -4.71721 1.54679 28.71297 -3.050 0.004887 **
## WindSpeed_KMI -0.20065 0.05343 143.30993 -3.755 0.000251 ***
## Temperature_KMI 0.02873 0.03580 78.77035 0.803 0.424623
## Cloudcoverage_KMI -0.17268 0.21176 142.63991 -0.815 0.416163
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Trj100 WS_KMI Tm_KMI
## Traject100m -0.632
## WindSpd_KMI 0.385 -0.266
## Temprtr_KMI -0.897 0.428 -0.593
## Cldcvrg_KMI -0.343 0.097 -0.336 0.322model2<-lmer(ForagingSpeed ~ Traject100m + WindSpeed_KMI + Cloudcoverage_KMI + (1|NestID) + (1|NestID:BaitID) + (1|Individualcode), na.action=na.omit, data=Flights2)
summary(model2)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: ForagingSpeed ~ Traject100m + WindSpeed_KMI + Cloudcoverage_KMI +
## (1 | NestID) + (1 | NestID:BaitID) + (1 | Individualcode)
## Data: Flights2
##
## REML criterion at convergence: 465.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.47342 -0.52496 0.04095 0.44893 2.92491
##
## Random effects:
## Groups Name Variance Std.Dev.
## Individualcode (Intercept) 0.1711 0.4136
## NestID:BaitID (Intercept) 0.9100 0.9539
## NestID (Intercept) 0.2484 0.4984
## Residual 0.3224 0.5678
## Number of obs: 191, groups: Individualcode, 74; NestID:BaitID, 52; NestID, 23
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 4.59898 0.32492 32.71481 14.154 1.67e-15 ***
## Traject100m -5.27920 1.37511 28.50382 -3.839 0.000632 ***
## WindSpeed_KMI -0.17545 0.04319 132.44886 -4.063 8.27e-05 ***
## Cloudcoverage_KMI -0.22488 0.20110 140.07082 -1.118 0.265372
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Trj100 WS_KMI
## Traject100m -0.620
## WindSpd_KMI -0.422 -0.017
## Cldcvrg_KMI -0.131 -0.047 -0.191model3<-lmer(ForagingSpeed ~ Traject100m + WindSpeed_KMI + (1|NestID) + (1|NestID:BaitID) + (1|Individualcode), na.action=na.omit, data=Flights2)
summary(model3)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: ForagingSpeed ~ Traject100m + WindSpeed_KMI + (1 | NestID) +
## (1 | NestID:BaitID) + (1 | Individualcode)
## Data: Flights2
##
## REML criterion at convergence: 557.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.56525 -0.54882 0.06406 0.46203 3.02412
##
## Random effects:
## Groups Name Variance Std.Dev.
## Individualcode (Intercept) 0.1638 0.4047
## NestID:BaitID (Intercept) 0.7251 0.8515
## NestID (Intercept) 0.3533 0.5944
## Residual 0.3176 0.5636
## Number of obs: 230, groups: Individualcode, 91; NestID:BaitID, 68; NestID, 34
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 4.6054 0.2843 63.4475 16.197 < 2e-16 ***
## Traject100m -5.1279 1.2857 41.4125 -3.989 0.000265 ***
## WindSpeed_KMI -0.1895 0.0394 169.3697 -4.810 3.31e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Trj100
## Traject100m -0.606
## WindSpd_KMI -0.517 -0.011library(cAIC4)## Loading required package: stats4## Loading required package: nlme##
## Attaching package: 'nlme'## The following object is masked from 'package:lme4':
##
## lmList## The following object is masked from 'package:dplyr':
##
## collapselibrary(MASS)##
## Attaching package: 'MASS'## The following object is masked from 'package:plotly':
##
## select## The following object is masked from 'package:dplyr':
##
## selectlibrary(AICcmodavg)##
## Attaching package: 'AICcmodavg'## The following object is masked from 'package:lme4':
##
## checkConvAIC(fullmodel)## [1] 152.8988AIC(model1)## [1] 487.598AIC(model2)## [1] 481.3893AIC(model3)## [1] 571.1946